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^ '. ABSTRACT 

, This paper investigates the sound production in a system made of a bore coupled with a reed 

■ valve. Extending previous work (Debut, 2004), the input impedance of the bore is projected on 
^ I the modes of the air column. The acoustic pressure is therefore calculated as the sum of modal 
Q ■ components. The airflow blown into the bore is modulated by reed motion, assuming the reed to 

•^1 be a single degree of freedom oscillator. Calculation of self-sustained oscillations controlled by 

time-varying mouth pressure and player's embouchure parameter is performed using ode solvers. 
^ [ Results emphasize the participation of the whole set of components in the mode locking process. 

, Another important feature is the mutual influence of reed and bore resonance during growing 

blowing pressure transients, oscillation threshold being altered by the reed natural frequency and 
the reed damping. Steady-state oscillations are also investigated and compared with results given 
^ ' by harmonic balance method and by digital sound synthesis. 
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(-^j ■ INTRODUCTION 

One of the main goals of musical acoustics is to provide models and investigation methods to 

■ improve our understanding of the behaviour of musical instruments, with application to the sound 
I synthesis and assistance of instrument makers. In the case of woodwind instruments, knowledge 

O ■ about sound production mechanisms has improved with models including the intrinsic nonlin- 

ear coupling with the excitation. Steady-state oscillations have been investigated in frequency 
domain [1 1 , 5], whereas time domain calculations were proposed since [Wi for characterizing 
transients and stability of oscillations. 
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The aim of this paper is to present a calculation of self oscillations in time domain based on the 
modal decomposition of the bore input impedance. Modal analysis is widely used in musical 
acoustics to analyse and reproduce the vibration of complex vibrating structures, but few appli- 
cations have been presented for self-sustained instruments, [1] for the study of the bowed string, 
and [9] for sound synthesis. The outline is as follows : the next section presents the model of a 
clarinet-like system and the formulation adopted to compute self-sustained oscillations. We then 
present the results of calculations, and discuss the behaviour of the different components of the 
pressure field, and the role of reed dynamics on transients. Finally this method is compared to 
sound synthesis algorithm and harmonic balance technique. 



MODEL AND IMPLEMENTATION 
Reed motion 

The elementary model, already used by many authors |8], relates the volume flow entering the 
mouthpiece u{t), and the pressure in the mouthpiece p{t) to the reed displacement y{t) by means 



of three equations characterizing the reed motion, the resonator and the flow of air through the 
reed channel. We assume for the characterization of the reed motion a lumped second-order 
mechanical oscillator with stiffness K per unit area, damping q,. and natural angular frequency ujr, 
driven by the pressure drop Pm - pit) accross the reed with an inward striking behaviour : 

Pm being the blowing pressure. This assumption, based on the fact that reed displacement occurs 
in the vertical direction mainly without torsion, has been discussed [2J, and appears to be valid for 
non-beating reed regimes. 



The airflow 

As noted by Hirschberg [7J, in the case of clarinet-like instruments, the control of the volume flow 
by the reed position is due to the existence of a turbulent jet. Indeed, a jet is supposed to form 
in the embouchure (pressure pjct) after the flow separation from the walls, at the end of the (very 
short) reed channel. Neglecting the velocity of air flow in the mouth compared to jet velocity t^jct, 
the Bernoulli theorem applied between the mouth and the reed channel leads to: 

P,n = Pjct + ^pwjct when p is the air density. (Eq. 2) 

Assuming a rectangular aperture of width W and height H + y{t) (H being the height without any 
pressure difference), the volume flow u accross the reed channel can be expressed as follows: 

uit) = WiH + y{t))^^^Prr,-piet{t). (Eq. 3) 

Since the cross section of the embouchure is large compared to the cross section of the reed 
channel, it can be supposed that all the kinetic energy of the jet is dissipated through turbulence 
with no pressure recovery (like in the case of a free jet). Therefore, the pressure in the jet is 
(assuming pressure continuity) the acoustic pressure p{t) imposed by the resonator response to 
the incoming volume flow u. 



Input impedance of the bore 

We consider a cylindrical bore (with length L) for the clarinet. The retained model - described 
extensively in [3] - is the wave equation inside the tube, with Neumann and Dirichlet boundary 
conditions at the input and output of the bore respectively, and a source taking into account the 
airflow blown inside the instrument: 

{9l^-{'^ + '^)^)Pi'^^x)^-juj§U{uj)S{x-Xs) yxe[0,L] andxs^O, 
da:P{uj, x)^0 for X = 0, (Eq- 4) 

P{uj, x) = for X = L. 

The pressure field p{x, t) is expanded onto the modes of the air column inside the bore: p{x, t) = 
Ht^i fn{x)pn{t) where the family {/n}„gN is a basis of orthogonal eigenmodes. In the case of a 
close/opened cylindrical pipe, fn{x) = cos (fc„a;), where fc„ = n-KjlL and n is an odd positive num- 
ber. Modal coordinates pnit) are calculated through the projection on each mode /„ of equation 
l |Eq. 4| in the time domain, the expansion onto the modes being truncated to the TV first modes: 

% + 2a„c%i + fc^cV(t) = ^^c$ for n £ [1, N\ and = ^. (Eq. 5) 

dt^ at L at b 

As explained in [3], we consider that the value a of damping associated to each mode may be 
different and is frequency-independent in the neightboorhood of each angular frequency w„ = fc„c. 
Damping is then noted ~ Yn/L where F„ is the value of the admittance at 
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Dimensionless model 

Considering as a reference tine pressure pm = KH required to close tine reed cinannel in non- 
oscillating case, we introduce here the dimensionless quantities for pressure in the mouth and in 
the mouthpiece and for volume flow: 

-P(^) - ZrU{t) 

pit)^ , Pn{t)^ , 7 = , u{t)^ . (Eq. 6) 

Pm Pm Pm Pm 



Likewise it is convenient to define x{t) = y{t)/H + 7 and C = ZcWH ^j2{ppMy^ , such that in 
the case of constant control parameters 7 and C, we have a complete model made of TV + 1 
second-order ode: 



x dx 2 / \ 2 ~/ 



^ + 2a„c^ + (^^-«VK(t) = ^§ fornG[l,iV], (Eq. 8) 



with u = ({1 + X - -f) sign (7 - p) v^|7 (Eq. 9) 

It is important to enhance that parameters 7 and (, controlled by the player while playing, are 
allowed to vary in time. This time dependency involves additional terms in the previous equations 
which allows the simulation of initial transients with a growing blowing pressure. The complete 
system is not written here but results are discussed in the next section. 



Beating reed 

Avanzini et al [2] investigated the effect of reed curling and beating against the lay of the mouth- 
piece, and proposed a model of variable stiffness. We include in our model the modelling of 
contact between lay and reed by considering an additional elastic force Fiay = Kigy{H + y{t)) 
when the tip gets in contact with the lay and completely closes the channel . Further investiga- 
tions would carry the case where the contact point is not the tip : as the reed curls against the 
mouthpiece, stiffness increases as the part of the reed which can freely move decreases, and 
then grows abruptly when the tip hits the lay. Another consequence of reed beating is that the 
volume flow vanishes when the reed closes the channel: the contact force prevents the reed from 
going beyond the lay in simulation, and so ensures that l+x-^ > 0, provided /^lay is high enough. 
This coefficient is here assumed to be an tunable parameter due to the lack of experimental in- 
vestigation in this domain. Typical value is about 100 times the stiffness of the non-beating reed. 



Numerical formulation 

The set of equations is processed using ode solvers. It constrains the formulation of system 
l |Eq. 7|Eq."9) to be rewritten as a set of 2(A^ + 1) first order ode, ( |Eq. 9^ being used in the source 
term of ( |Eq. 8| . As multiple time scales are involved in the problem (duration of blowing pressure 
transient, bore resonance frequencies, reed resonance frequency and reed-lay contact duration), 
we used several ode solvers designed for stiff problems, in particular odelSs from the Matlab ode 
Suite and isode from odepack, based on the Backward Differentiation Formula (bfd) methods. 
They lead to similar results in typical cases, for example a difference for the oscillation fundamental 
frequency calculated of only 2 cents for a bore length equal to 0.5m (oscillation fundamental 
frequency : /o = 170Hz). 



RESULTS 

Simulation of attack transients 

Using the possibility to control the mouth pressure, simulations of attack transients with increasing 
blowing pressure are carried out. The time evolution of 7 is as follows: 

{0 for t < 0, 

7o/2(l -cos(7rt/T)) fort > Oand t < r, (Eq. 10) 

70 for t > T. 
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Figure 1 : Pressure p{t) and components {pn{t)}n=i,8 for two attack transient simulations: on tine 
left side, instantly growing mouth pressure, and on the right, slowly growing blowing pressure 
(10ms) with qr = 0.4, cjr ^ 2n x 1500Hz for both simulations. Waveforms in steady-state for each 
signal are also displayed. 



Figure [H shows two simulations : on the left picture, a quasi instantaneous increase of mouth 
pressure (r = 1/44100 s) and in the right one a slower increase (r = 10ms). The difference 
between the results have to be highlighted. On the right figure, blowing pressure slowly varying 
from to a value high enough to make the oscillation start: component pi is the first to appear, 
stabilizing its amplitude by the nonlinear phenomena in coupling. Then this saturation produces 
the apparition of the others components (since t = 0.05s). With an abruptly established mouth 
pressure (left figure), the whole set of components are excited, but there is a differentiation be- 
tween decaying components and increasing ones. These latter are called master components, 
and are involved in self-sustained oscillations generation. On the other hand, linear stability anal- 
ysis allows to conclude to the stability of decaying component for such blowing pressure, but 
nonlinear coupling implies energy transfer between components and the initially decaying ones 
finally converge to a non-zero oscillation: there are the slave components [3]. 

Another interesting feature is the fact that higher components seems to have their amplitude 
modulated by fundamental oscillation (see waveforms in figure[T), showing the difference between 
time domain simulation involving modal decomposition and a simple signal decomposition based 
on Fourier analysis: each first order resonator is excited by the input volume flow variation, that is 
by a rich spectrum. This results to the fact that each component is not a monochromatic signal. 



Comparison with l-iarmonic balance method 

Figure|2]displays the result of the comparison between two time domain ode solvers ( odel5s from 
the Matlab ode suite and isode from odepack) and harmbal (a software for computation of pe- 
riodic solutions by the harmonic-balance method, [4]) for a periodic regime. Differences between 
HARMBAL and time-domain solvers appear to be of the same order as differences between the 
two time-domain solvers. However, unlike time-domain methods, the harmonic balance method 
does not capture transients, therefore only periodic states can be compared. 
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Figure 2: Pressure p power spec- 
trum calculated with different solvers 
showing good agreement (7 = 0.45, 
C = 0.4). Solvers considered include 
odel5s from the Matlab ode suite, 
isode from ODEPACK and Harmbal 
from [4]. 



Comparison with sound synthesis 

The time domain scheme used for digital sound synthesis consists in obtaining digital filters rep- 
resenting both the reed displacement fEqTTj and the input impedance. The reed displacement 
is discretized using centered finite difference schemes for the first and second order derivatives. 
This leads to a non-instantaneous response of the reed filter (i.e y[n] does not depend on p[n] but 
on p[n - 1]). The impedance is modeled by a digital filter expressing, at sample n, the pressure 
p[n] as: p[n] = bcou[n] + V[n], where V[n] includes past (and known) values of p and u. For 
the purpose of this paper, this digital filter is computed from the inverse Fourier transform of the 
reflexion coefficient associated to the modal impedance. These discretizations allow to solve ex- 
plicitely the nonlinear coupling, whatever the bore geometry, as it has been shown in [6]. Since the 
modal impedance comes from a truncated serie of modes and does not take into account radiation 
losses or tonehole network, it is worth noting that at high frequencies, the associated reflection 
coefficient tends to 1 and the admittance tends to infinity. This behavior generates instabilities in 
the acoustic flow during attack transients and we are working on this problem. 
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Figure 3: Pressure p signal resulting 
from sound synthesis (top picture) and 
simulation (bottom picture) for a in- 
creasing blowing pressure (70 = 0.45 
and r = 10ms). 



CONCLUSION 

The approach presented in this paper relies on a time domain formulation, leading to ab ini- 
tio calculations. Therefore, different stable regimes are computable depending obviously on initial 
conditions but also on the possibly time-varying values of control parameters. Moreover the modal 
decomposition of the dynamics allows an explicit writing of the problem using meaningful param- 
eters in musical acoustics such as natural resonance frequencies and dampings. Since these 
quantities are accessible to measurements, our simulations can be reasonably easily compared 
to experiments. This is not always the case with time domain methods using impulse responses 
or reflection functions for example [10]. 

Another advantage of the time-domain modal approach is that the resulting formulation (a system 
of first-order odes) is the same as the one used in dynamical system studies. Benefits are ob- 
vious to analyze behaviors observed in simulations (for example the analysis of attack transients 
in terms of master/slave components in section Results of this paper). A more practical conse- 
quence is that many reliable solvers are available for first-order ode systems. This allows to test 
many routines without additional development, in order to check the validity of a simulation result, 
or in order to look for the best solver for a particular problem. For example, this appeared to be 
decisive when considering reed beating regimes, very demanding numerically, which required a 
routine specifically designed for stiff problems. 

The computational cost may appear high compared to other implementations. This is not a limi- 
tation of the time domain modal formulation since other similar approaches are used for real-time 
sound synthesis [9J, but is mainly due to the type of chosen solver (variable order, multi-step 
scheme). Our approach is thus definitely not intended for sound synthesis. On the other hand, 
our simulations are expected to be more precise and aimed at investigating model dynamics. 
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